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We derive an effective equation of motion for the orientational dynamics of a neutrally buoyant 
spheroid suspended in a simple shear flow, valid for arbitrary particle aspect ratios and to linear 
order in the shear Reynolds number. We show how inertial effects lift the degeneracy of the Jeffery 
orbits and determine the stabilities of the log-rolling and tumbling orbits at infinitesimal shear 
Reynolds numbers. For prolate spheroids we find stable tumbling in the shear plane, log-rolling is 
unstable. For oblate particles, by contrast, log-rolling is stable and tumbling is unstable provided 
that the aspect ratio is larger than a critical value. When the aspect ratio is smaller than this value 
tumbling turns stable, and an unstable limit cycle is born. 

PACS numbers: 83.10.Pp,47.15.G-,47.55.Kf,47.10.-g 


I. INTRODUCTION 


In this article we describe the effect of weak inertia upon the orientational dynamics of a neutrally buoyant 
spheroid in a simple shear flow using perturbation theory. In the absence of inertial effects the rotation of a 
neutrally buoyant spheroid in a simple shear was determined by Jeffery who found that there are infinitely 
many degenerate periodic orbitJi^, the so-called ‘Jeffery orbits’. In this limit the initial orientation determines 
in which way the particle rotates. Fluid and particle inertia lift this degeneracy, but little is known about 
how this comes about. A notable exception is the work by Subramanian and Koch who have solved the 
problem for rod-shaped particles in the slender-body approximatiorP. We discuss other theoretical results 
below in Section m 

The question is currently of great interest: several recent papers have reported results of direct numerical 
simulations (DNS) of the problem, using ‘lattice Boltzmann’ methodP^^. These studies reveal that fluid 
and particle inertia affect the orientational dynamics of a neutrally buoyant spheroid in a simple shear 
in intricate ways. The DNS are performed at moderate and large shear Reynolds numbers, defined as 
RBs = sc? jv where a is the largest particle dimension, s is the shear strength and v the kinematic viscosity 
of the suspending fluid. DNS at very small Reynolds numbers are difficult to perform. But this limit (Reg 
of order unity and smaller) is of particular interest. There is a long-standing question whether or not a 
nearly spherical prolate spheroid exhibits stable ‘log-rolling’ in this limit, so that its symmetry axis aligns 
with the vorticity axis. It was first suggested by Saffman that this is the casP, in an attempt to explain 
Jeffery’s hypothesiP that spheroids rotate in orbits that minimise energy dissipation. But stable log-rolling 
of prolate spheroids has not been found in DNS, and it has been suggested that higher Reg-corrections may 
explain this discrepancjP. The small-Reg limit is of interest also because it provides stringent tests for DNS. 
These reasons motivated us to derive an equation of motion that takes into account the effect of weak fluid 
and particle inertia. Our main result is an approximate dynamical equation for the rotation of a neutrally 
buoyant spheroid suspended in a simple shear flow, valid for arbitrary aspect ratios and to first order in Reg 
(Eq. in Section IV). In the slender-body limit this equation is of the same form as the one derived in 
Ref.[^ In the completely inertia-free case our results reduce to Jeffery’s equatioiP. We find that corrections 
to this limit arise from both particle inertia (centrifugal and gyroscopic forces), as well as from fluid inertia 
(modifying the hydrodynamic torque on the particle). The particle-inertia corrections we report here are 
consistent with earlier numerical and analytical resultd^. 

Fluid-inertia corrections are taken into account to first order in Reg using a reciprocal theoren^. Our 
approach is similar to the one adopted in Ref. [2] in the slender-body limit, but our equation of motion is 
valid for spheroids with arbitrary aspect ratios. By linear stability analysis we determine the stabilities 
of the periodic orbits of this equation at infinitesimal Reg as a function of the particle aspect ratio. The 
stability calculation details how the degeneracy of the Jeffery orbits for a neutrally buoyant spheroid in a 
simple shear is lifted by weak inertia. 

We find that the log-rolling orbit is unstable for prolate particles. This explains why stable log-rolling is 
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FIG. 1. Color online. Spheroid rotating in a simple shear (schematic). Shows the Cartesian coordinate system ej, 
j — 1,..., 3. Vorticity points in the negative ea-direction. The flow-shear plane is spanned by ei and ea. The unit 
vector n points along the symmetry axis of the spheroid. Its polar angle is denoted by 9, the azimuthal angle by ip. 
The major axis length of the spheroid is denoted by a, the minor length by b. For prolate particles the aspect ratio 
is defined as X = a/b, and for oblate particles as A = 6/a. 


not observed in the smallest shear Reynolds numbers accessible in the simulations. Moreover we 

find that tumbling in the flow-shear plane is stable for prolate particles. As the aspect ratio tends to unity 
there is a bifurcation: for nearly spherical oblate particles log-rolling is stable and tumbling in the flow-shear 
plane is unstable. There is a second bifurcation for oblate particles. At a critical aspect ratio Ac « 1/7.3 
tumbling becomes stable and an unstable limit cycle is born. This means that the behaviour of a very flat 
disk depends on its initial orientation for A < Ac. We discuss how the shape of the limit cycle changes as 
the aspect ratio tends to zero. 

The remainder of this art icle is organised as follows. In Sectionj^we give an overview over the background 
of the problem. Section III summarises the method employed in this article, based on a reciprocal theorerrfi^. 
We demonstrate how to calculate the effect of particle and fluid inertia to first order, and how we use the 
symmetries of the problem to make it tractable. Section [TV] summarises our results: the equation of motion 


and its stability analysis. We discuss the results in Section [Vj and conclude with Section [VT 


A brief account of the main results described in this article was given in Ref. [TI] Here we describe the 
complete derivation. We also present additional results and discussion that could not be included in the 
shorter format: we quote precise asymptotic formulae for small and large aspect ratios, as well as for aspect 
ratios close to unity. We also characterise the limit cycle that arises for A < Ac, and compute its linear 
stability. 


II. BACKGROUND 

The question of describing the rotation of a neutrally buoyant particle in a simple shear flow has a long 
historv Jeffery derived an expression for the torque on an ellipsoidal (tri-axial) particle neglecting inertial 
effect^. To obtain an equation of motion for small particles he assumed that the dynamics is overdamped, 
that the particle rotates so as to instantaneously achieve zero torque. This gives rise to Jeffery’s equation that 
is commonly quoted for the special case of spheroidal (axisymmetric) particles. From this equation it follows 
that spheroids suspended in a simple shear tumble, they stay aligned with the flow direction for some time 
and then switch orientation by 180 degrees. The dynamics is degenerate in that there are infinitely many 
different periodic orbits, the so-called ‘Jeffery orbits’. The initial orientation determines which particular 
orbit is selected. The goal of Jeffery’s calculation was to compute the viscosity of a dilute suspension of 
spheroids, and Jeffery hypothesised that the particles select orbits that minimise energy dissipation. 

SaffmarP pointed out that inertial effects lift the degeneracy of the Jeffery orbits, and he described the 
orientational dynamics of a nearly spherical particle in a simple shear taking into account fluid inertia. For 
prolate particles he concluded that log-rolling is stable, that tumbling in the shear plane is unstable, and 
that the stabilities are reversed for oblate particles. These results are stated in terms of an effective drift for 
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the particle orientation (towards the vorticity axis for prolate particles). This conclusion supports Jeffery’s 
hypothesis. Saffman did not take into account particle inertia. His method of calculation rests on a joint 
expansion in small eccentricity and RCg. 

Harper & Chan^i^ addressed the problem in a different way, modeling the dynamics of a rod in a simple 
shear in terms a dumb-bell, that is two spheres connected by an invisible rigid rod. The spheres are subject 
to Stokes drag and hydrodynamic lift forceJ^^. This approximation neglects hydrodynamic interactions 
between the two spheres, as well as the unsteady term in the Navier-Stokes equations. Harper & Chang 
arrive at the opposite conclusion, namely that log-rolling is unstable. Since their result pertains to the 
slender-body limit the question is how the stability of the log-rolling orbit depends on the particle aspect 
ratio. 

It was subsequently shown by Hinch & LeaP^ how weak rotational diffusion breaks the degeneracy of the 
Jeffery orbits, and their results form the basis for a large part of the work during the last decades on the 
rheology of dilute suspensions, see Refs. [T3] and M for reviews. 

Recently there has been a surge of interest in determining the effect of weak inertia upon a spheroid 
tumbling in a simple shear flow in the absence of rotational diffusion. Subramanian & KoclP derived an 
effective equation of motion for a neutrally buoyant rod in the slender-body limit to first order in fluid and 
particle inertia. Their calculation uses a reciprocal theorenP^ and takes into account the unsteady term in 
the Navier Stokes equation as well as particle inertia. The authors arrive at qualitatively the same conclusion 
as Harper & Chang, namely that the orientation of the rod eventually drifts towards the flow-shear plane. 

In a second paper Subramanian & KoclP^ repeated Saffman’s calculation for a neutrally buoyant nearly 
spherical particle. They used a different method, similar to the one used in Ref. [21 and come to the same 
conclusion as Saffman, that log-rolling is stable for nearly-spherical prolate particles. 

Recent have explored the stability of log-rolling and tumbling orbits, mostly at moderate and 

large Reynolds numbers, and only for a small number of aspect ratios. The simulations show unstable 
log-rolling for prolate particles at the smallest Reynolds numbers studied. We note that You, Phan-Thien 
& TanneP^ misquote Saffman when they describe their numerical results on the rotation of a spheroid in a 
Couette flow at Reynolds numbers of the order of 10 and larger. In the introduction of Ref. [TH]it is implied 
that Saffman’s theorjPl predicts that nearly spherical prolate particles tend to the flow-shear plane. 


III. METHOD 


In this section we give a brief but complete summary of our calculation. The most technical details and 
tabulations are deferred to appendices. We start with notation and the relevant dimensionless parameters 
determining the physics. Then we giv e the gove rning equations and explain how to express the hydrodynamic 
torque through a reciprocal theorenpHSMEo] penally we explain the perturbation scheme and list the 
symmetries that severely constrain the form of the solution. 


A. Notation 

The calculations described in this paper involve vectors and tensors in three spatial dimensions. We employ 
index notation with the implicit summation convention for repated indices, and we use the Kronecker 
and Levi-Civita (e^fc) tensors. 


B. Units and dimensionless numbers 


The physics of the problem is governed by three dimensionless numbers: the shear Reynolds number Res 
(measuring fluid inertia), the Stokes number St (measuring particle inertia) and the particle aspect ratio A. 

We work with dimensionless variables. The length scale is given by the particle major axis a. The velocity 
scale is taken to be sa where s is the shear rate. The explicit time dependence of the flow (the time scale 
for the unsteady fluid inertia) scales as ~ 1/s since it is determined by the particle angular velocity because 
to lowest order the unsteadiness arises from the particle motion. The corresponding scale for pressure is ^s, 
and force and torque are measured in units of and /rso^, respectively. 
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From these scales the dimensionless parameters are formed. As mentioned in the Introduction, the shear 
Reynolds number is defined as 


Res 


so? Pi 


( 1 ) 


Here pf is the density and p = piV is the dynamic viscosity of the surrounding fluid {v is the kinematic 
viscosity). 

The Stokes number, measuring the particle inertia, is given by the ratio of the typical rate of change of 
angular momentum and the typical torque: 


St=^ 


^Res. 

Pf 


( 2 ) 


Here Pp is the particle density. For a neutrally buoyant particle, Pp = pf, we have St = Re^. 

We define the particle aspect ratio A as the ratio between the length along the symmetry axis and the 
length transverse to the symmetry axis (Fig. [^. That is A > 1 denotes prolate particles while A < 1 denotes 
oblate particles. Because we measure length in units of the major particle axis a, the aspect ratio A of a 
prolate particle is a/&, while the aspect ratio of an oblate particle is 6/a, where b denotes the length of the 
minor axis of the particle (see Fig. [^. 


C. Equations of motion 


Let rii denote the components of the unit vector n pointing in the direction of the particle symmetry 
axis (Fig. [^, and uJt the components of the angular velocity of the particle. Newton’s second law for the 
orientational degrees of freedom for an axisymmetric particle reads: 


^2 — ^ijk^j‘^k 5 


St 


lijUJj 




= T,,. 


( 3 ) 


Dots denote time derivatives, and lij are the elements of the moment-of-inertia tensor of the particle, and 
Tj is the torque exerted on the particle. The moment-of-inertia tensor of an axisymmetric particle with axis 
of symmetry n is on the form 


lij = A^Tiirij + {Sij — riirij). 


( 4 ) 


where and correspond to the moments-of-inertia around and transverse to the symmetry axis. Using 
the dimensionless variables introduced in Section HIB we have for a prolate spheroid (A > 1) 


A^ = 


Stt 1 

IsF’ 


B^ = 


4tt 1 
15 W 


1 + 


A2 


( 5 ) 


and for an oblate spheroid (A < 1) 


^ "15'^’ 


We rewrite the equation of motion ^ as 


47r 


^^ = iA(l + A^). 


— hj Ijk^k + Tj — 


A^-B^ 1 


B^ 


( 6 ) 


( 7 ) 


In the final step we used the definition Q of lij and the equation of motion ^ for hi. 

In this paper, the Ti are the elements of the hydrodynamic torque exerted on the particle by the fluid. In 


Section HID we formulate the hydrodynamic torque to O(Res) via the reciprocal theorem. In Section HIE 
we perturbatively compute the resulting angular velocity to order O(Res) and 0(St). 
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D. Calculation of the hydrodynamic torque to order Res 


The straightforward approach to determine the torque on a particle in a fluid is to solve Navier-Stokes 
equations for the velocity and pressure fields, then compute the stress tensor and Anally integrate the stress 
tensor over the surface of the particle. The reciprocal theorenPEHIIIl offers an alternative, and often more 
convenient, route to the hydrodynamic forces. In particular, we may avoid solving for the complete flow 
field. In this Section we specify the Navier-Stokes problem we need to solve, and explain how we use the 
reciprocal theorem to simplify the calculations. 

Navier-Stokes problem for the disturbance flow. We consider a particle with boundary S immersed in an 
linear ambient flow Throughout this paper we express the ambient flow as 


— A°°r. 


= Ei 


ik3^‘‘k 'j 


S°°r- 

^13 '3 


or equivalently with Sikj^'^ = 0°° 


( 8 ) 


uT = Of°r, + Sf°r,. 


Here and are the symmetric and antisymmetric parts of the flow gradient, given by 


ooo _ ^ A^\ _4'^^ 

^ij — 2 ’ ^^3 ~ 2 


In dimensionless variables (Section III B[) the Navier-Stokes equations read 


Reg ~\~ UjdjUi^ — dip H- djdjUi , diUi 


= 0 . 


( 9 ) 


( 10 ) 


( 11 ) 


Note that the unsteady and convective inertia terms come with the same prefactor in this problem. This 
happens because the timescale of the particle motion is the same as the timescale of the flow. The boundary 
condition is no-slip on the surface of the particle 


Ui=eijkOJjrk toi rGS , Ui=uf° as|r|—>- 00 . (12) 

We introduce the disturbance held {ul,p') from the particle 

Mi = -k u', p = -k p'. (13) 

If we assume that {u°°,p°°) satisfies the Navier-Stokes equations we have the disturbance problem 

Re^ {dtUi ufdju'^ u'^djuf -\- u'pju^) = —dip' djdju [, (14) 

and the boundary conditions is expressed in the slip angular velocity Hi = — uji as 

Mi — Eijk^jrk rj.j V (z: Sj 

m( = 0, |r| —>■ oo . (15) 

Finally, when applying the reciprocal theorem we shall use that, by definition, the divergence of the stress 
tensor satisfies the following equalities: 

dja'ij = -dip' + djdju'i 

= Reg {dtu[ -\- u^dju[ u'jdju^ u'pju'f) 

= Resfi{u'). (16) 


The Stokes solution. T his paper concerns a spheroidal particle suspended in a linear flow. We thus need 
explicit solutions to Eq. (14) at Re^ = 0 in this geometry. We use a finite multipole expansiorP^III] (see 
Appendix [A| . In our notation they read 


— Qij,k^jkl [(-^ T ^3 riiTlyyif^ 4” ^ Ei.fym’^mS.^QriQ^ 

X [(A B -j- C 5";^ — C {Sjlmnkrim + Eklmnjnm) ; 


( 17 ) 
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where 


— {’njUk ^Sjk)(ninm 

^jklm ~ kljSklklm “t“ "klji^SjlTlrn “t“ TljSkrnk^l klk^jm’k^l '^'k^j’klkkilklm: 

^jklm ~ ^jk^lm “t“ Sjl6km “t“ ^kl^jm 

-\- Sjf^TllTlYfi -\- Sl'^TljTlf^ Tlj5]^lTlrui 

klj^km'kil 'k^k^jm^l “ 1 “ 'kljklk'l^lklm- 

Here A^, B^, C^, and a are known constants that depend on the particle aspect ratio A. The 

exact definition of the spheroidal multipoles Q and the values of all constants are given in Appendix]^ see 
in particular Table |III| 

The reciprocal theorem. This theorenPHnUlIini relates integrals of the velocity and stress fields of two 
incompressible and Newtonian fluids. The idea is the following. Let one set of fields represent the actual 
problem of interest, the primary problem. Then choose the second set of fields to be an auxiliary problem 
with known solution, such that an integral in the theorem relates to hydrodynamic torque of the primary 
problem. Provided that all integrals in the theorem converge and can be evaluated, we can solve the resulting 
equations for the hydrodynamic torque. 

The reciprocal theorem for the two sets {ui,dij) and {u[,a[j) can be stated as 


dib-i 


IV 


dVu'^djdij = 


dF'ui + 


IV 




(18) 


Here dFi = aij^jdS is the differential force from the fluid on the surface element with normal vector ^jdS. 
The volume integrals are to be taken over the entire fluid volume outside the particle, and the surface 
integrals over all surfaces bounding the fluid volume, with surface normals pointing out of the fluid volume. 

In the following we apply the reciprocal theorem to the calculation of the hydrodynamic torque on a 
particle. 

Calculation of the torque. We choose the auxiliary problem {ui,dij) to be the Stokes flow around an 
identical particle rotating with an angular velocity uji in an otherwise quiescent fluid. Its solution is given 
by Eq. @ with u°° = 0. The primary problem is the disturbance problem defined in Eq. ( |14[ ). Inserting 
the boundary conditions into the reciprocal theorem yields 


[ dFi{s^jkiuJj - - S^Tj) = [ dF'cijk^jrk + Re^ [ dVuifi{u'). (19) 

Js Js Jv 

We also used that djdij = 0. This equality holds because Ui is a Stokes flow. Both primary and auxiliary 
velocity fields vanish as |r| —>■ oo, therefore both integration surfaces are only the particle surface S. Note 
that the surface integrals are to be taken with surface normals out of the fluid domain, so that dFi is the 
differential force exerted on the particle by the fluid. In the integrals we identify the hydrodynamic torque 
on the particle, it is given by 


Ti = / dFiS 


ijk'^k • 


( 20 ) 


It follows: 


{ujj - )f, - / dF,S°°rj = Wj (T, - Tf) + Re, / dVu^Mu'). 

Js Jv 

The auxiliary torque Tj together with the surface integral add up to the Jeffery torqu^ : 

= cj [{A^UjUk + B^fSjk - UjUk)) (fl^ - Wfc) + C’^'SjkmnkniSf^i] . 

The constant cj is given in Table [Hi] in Appendix [A| The contribution 

P“ = J dSa°ff,ieijkCjjrk 


( 21 ) 


( 22 ) 


(23) 
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evaluates to zero for any linear flow It follows that Eq. (21) becomes 


ojjTj = — Res 


dVuif^iu'). 


Since Ui is linear in Cjj, this variable can be eliminated. We finally obtain: 


Tj = Tf - Res 


dVUijfi{u') , 


'v 


where 


(24) 


(25) 


U^p = -Qfj,k£oki [{A^niUp + B^{5ip - niUp))] (26) 

+ \p i^jpmnknm + SkpirinjUm)] ■ 

Thus far we have made no approximation, and Eq. ( |25[ ) is exact, the difficulty lies in evaluating the Navier- 
Stokes disturbance flow u'. This is a complicated non-linear problem since Uij and fi{u') all depend on 
the direction n and upon the angular velocity u) of the particle. The flow equations thus couple non-linearly 
to the rigid body equations of motion for the particle. In the following we solve this system of equations in 
perturbation theory valid to first order in St and Res. 


E. Perturbative calculation of the particle angular velocity 


In this section we determine the angular velocity ui of the particle to lowest order in St and Res, assuming 
that both St and Res are small, so that Reg St is negligible. We recall the equation of motion ([7|) for the 


particle orientation, and insert the expression for the hydrodynamic torque obtained in Section HID 

‘^i — ^ijk^j'^k 7 

o • -B^ 

StWi = -St— —p — Ei 

Now we expand the angular velocity as 


Next, we insert these expansions into the equation of motion ([^ and collect terms of equal order in St and 
Res: 


-ResI-/ l^dVUkjfkiu'). 

(27) 

^ -1- o(Res, St). 

(28) 


= - 


(0) 


- B^ 


(0) (0) r—l I aR I taR/s: \\ (St) 

EijkiA}) nkUiUJi - (A rijUk + B [djk - njUkjjujl , 


,(o) 


(St) 


\R, _I jjR/x . _W , ,(R® 


d = c^[A^njnk + B^{5jk-nink))^AJi‘+ / dVUkjfkiu'). 


(29) 


In the last term it is understood that the volume integral need only be evaluated to 0(1), so that we may 
use the Stokes flow solutions for it'. The first equation gives the Jeffery angular velocity uif'\ 


^ + ~^^ikmnkniS^i 


The dynamics of n is to lowest order given by 

fi^ — ^ Rq^ipq^p Rq ^R i^ip^P RiRpRq^pq') • 

From Table Ell in Appendix we infer that for both prolate and oblate spheroids 


(30) 


(31) 


_ A2 - 1 

' b ^~ ~ pTT ■ 


(32) 
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This shows that Eq. (31) is Jeffery’s equatiorfi^ for the orientational dynamics of a spheroid in a simple shear. 

may be inverted to 


The two remaining equations in 

(St) 1 
UJl = — 


AR 


Tl-iTlj RiRj') 


•,( 0 ) 


OJ. 


(Re.) 


1 / 1 




br^ y 

W 

1 


T (0) (0) 




BR' 


dVUkjfk{u')- 


(33) 


Eq. (28 1 together with Eqs. (301 and (331 yield the effective angular velocity under the effect of weak particle 
and fluid inertia. Erom the equation of motion ^ we define the effective vector field 


■hi = e^jkiOj'nk = + Sthf*^ 




(34) 


This vector field describes the time evolution of n. The first term is the Jeffery vector field (31). The two 
new terms represent the effects of particle inertia and fluid inertia. The terms due to particle inertia are 
straightforward to evaluate directly, but the volume integral in Eq. (33) is very tedious to evaluate. To make 
the calculation feasible we exploit the symmetries of the problem. 


F. Symmetries of the effective equation of motion 


Both correction terms in Eq. (341 are quadratic in the ambient flow gradient tensor 
they are on the form 


In other words. 


(35) 


where the tensorial coefficients are composed of the remaining available tensor quantities: rip and 5pq 

[sijk is already used in 0“)- We make an exhaustive enumeration of all possible combinations, and then 
use the symmetries listed in Table |T] to remove or combine items in the list. For example we start by letting 


^( 1 ) 

^ijklm 



T'P1^P2P3^P4:P5 


V 2 R‘Pl'R‘P2'R'P3^PiP5 


I (P) 

T Vs R-pi^p2^P3^P4^P5 


(36) 


where the sum is over all 5! = 120 permutations P of {i,j,k,l,m), and are unique coefficients for each 
term. We include only odd powers of rii as any even terms would break the particle inversion symmetry. We 
then insert this enumeration into the first term of Eq. (351, and contract and apply the first three symmetries 
in Table [T] until we reach a list of unique candidate terms. In this case the only two unique terms turn out 
to be OijOjkRk and ■nirijOjkOkini. Finally, we use the fact that the equation of motion may not change the 
magnitude of the unit vector n. This constraint forces the coefficients of the two unique terms to be the 
same magnitude but opposite sign. Upon renaming the coefficients ±ai we get the first term in Eq. (37). 
The other terms are derived similarly by inserting (36) into the other terms in Eq. (35). The result contains 
only six independent terms: 


hi = ai {6,j-n,'nj) O'^O'^ni 
+ 02 {dij-niUj) 

+ as {Sij-ni'Tij) 

+ 04 {Sij-riinj) 3^3^11.1 
+ as{np3'^qnq)0'^nj 


+ ais{jlp3pqPlq) {Sij TliPlj') 3jp.Tllz . 


(37) 


Here the scalar functions oi,..., 05 are linear in St and Re^, and depend on the aspect ratio 
linear (and unknown) way. These coefficients are determined by evaluating the vector field in 
six independent directions of n and solving the resulting system of linear equations for Oi,..., 


A in a non- 
Eq. ([34|) for 


06. 
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TABLE I. Symmetries constraining the form of the effective equation of motion, Eq. (34|. 


s!r = o 

Incompressible flow 

000 QOO 

Oij ‘^ji 

S°° symmetric 

o°° = -0°° 

0°° anti-symmetric 

riihi = 0 

Dynamics preserves magnitude 

Ui —Ui hi —hi 

Particle inversion symmetry 


In the particular case of a simple shear flow, we have explicitly (see Fig. [^for the geometry) 

i iSaS,2-6M . ^ iS^lS,2 + 6,2S,l). (38) 

We observe that for the simple shear 0^0^ = S^S’^, and Sf^O^ = -Of°S^. Then the form of the 
equation of motion simplifies to 

ni — [ 5 \{lflpSpgTiq) {^Sij TliTlj^ 

+ P2inpS^gng)0^nj 
+ Pa {S^J - n^Tij) Of^S^iUi 

+ Pi {6,j - n^Tij) S^S^ni, (39) 


with 


Pi = CtQ, 

P 2 = , 

Pa = ct3 — o :2 , 

Pi = ai — ai . (40) 


Thus, for the case of the simple shear flow it suffices to evaluate the effective vector field in Eq. (34), in 
particular the volume integral in Eq. (33), with four independent values of n in order to solve for the 
unknown scalar coefficients Pi^... ^ Pi. 


G. Evaluation of the volume integral in Eq. (|33[) 


The volume integral in Eq. (33) contains four distinct terms: represents unsteady fluid inertia, and 

the three terms u^dju[ + UjdjU^ + u'^ dju’^ represent convective fluid inertia. We compute these four terms 
using the explicit Stokes-flow solutions (17). While the Stokes flow has no explicit time dependence, both 
particle direction n and angular velocity w do. Thus each occurence of and uik has to be differentiated 
to compute the contribution due to unsteady fluid inertia. The differentiation and tensor contractions are 
implemented by a custom set of pattern matching rules in Mathematica®. The calculation is both long and 
error prone. We have therefore automated every possible step, including solving the Stokes-flow equations. 

We demonstrate the remainder of the procedure by a small example. Consider the contribution in the 
ea-direction of Eq. (33) due to unsteady fluid inertia: 


—5ia — 


ar 


niUj 


BR 


ipij 


Res 

St 


dVUkjdtu^. 


(41) 


We first perform the time derivatives on 0 in the manner explained above. Then we insert the components 
of n, and the explicit form of the shear flow ( [M| ). At this point we can explicitly perform the sum over all 
repeated indices. The result in this example consists of 858 terms, after collecting terms with same spatial 
dependence. The terms have a prefactor that stem from the Stokes-flow coefficients (see Appendix 0, 












10 


and a spatial dependence coming from and the spheroidal integrals and (see Appendix]^. For 
n = [1/2, •\/3/2,0] a typical term looks like this: 

1575a2C^(A‘5 - SC^){2B^ + 

16 (B«) 2 c^ 

We note that the only spatial dependence on the azimuthal angle around the symmetry axis of the body 
comes from factors r^. We introduce a rotated coordinate system in which = Rjir'j-, such that is 
along the particle symmetry axis (see Appendix . This change of basis enables integration of one spatial 
coordinate. 

After this operation 260 terms still remain which we program Mathematica® to express in spheroidal 
coordinates (Appendix]^ and integrate over the remaining two spatial coordinates. 

As a consistency check we have also evaluated the volume integral numerically over all three spatial 
dimensions by converting to spheroidal coordinates and choosing a specific value of A. For extreme values of 
A the numerics are difficult, nevertheless they serve as a check for a wide range of aspect ratios (see markers 
in Fig.[^. 


IV. RESULTS 

A. Effective equation of motion 

We parametrize the vector n in a spherical coordinate system (0, (p) with 9 the polar angle and p the 
azimuthal angle (Fig. [^: 


ni=sin6*cosvj 712 = sin 0 sin , 77,3 = cos 6*. 


In these coordinates Eq. (39) is expressed as 


ip[9, (^) = ^ (A cos 2p}— 1) + \j3\ sin^ 9 sin y sin 2p (/32 sin^ 9 + ^ , (42a) 

2 8 4 

9{9, (p) = A sin 9 cos 9 sin p cos p + ^ sin 0 cos 9 (/3i sin^ 9 sin^ 2p + cos 2p + . (42b) 

We compute the contributions to /3 q from three sources: particle inertia, unsteady fluid inertia and convective 
fluid inertia. Although the result is only valid for neutrally buoyant particles (Re^ = St), it is interesting to 
consider the contributions separately: 

y = Stpy + ResPy + ResPy (43) 

The contribution from particle inertia is straightforward to compute and can be expressed in closed form as 

(P) ^ 2B\Cy 

[BRfc^ ’ 

(P)_ - 2B^) 


P^2' = - 


pr = 


[Byc^ 


{BRfy 


(P) _ {A^ -B^){By +B^{Cy 


PI' = - 




(44) 


The coefficients on the r.h.s. of these equations are tabulated for both prolate and oblate spheroids in 
Table III in Appendix]^ The coefficients in Eq. (441 are shown as dotted lines in Fig. 

The expressions for the contributions from fluid inertia are very lengthy and not particularly instructive. 
We therefore present the full result graphically as function of aspect ratio A in Fig. In addition we give 
the asymptotic behavior of all contributions to Pa in Table jll] in three limiting cases: thin oblate particles 
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TABLE II. Asymptotic results for Pa- Contributions from particle inertia, unsteady fluid inertia, and convective 
fluid inertia are shown separately. Factors of Res and St are omitted. 


Thin oblate particles (A —>■ 0) 



Total 

Unsteady 

Convective 

Particle 

Pi 

11 

30 

1 

5 

1 

6 

0 

P 2 

1 

10 

1 

20 

3 

20 

0 

P3 

1 

5 

3 

20 

1 

20 

0 

P4 

1 

3 

3 

20 

11 

60 

0 


Nearly spherical particles (|e| <C 1) 



Total 

Unsteady 

Convective 

Particle 

Pi 

137e^ 

0 

163t'^ 

2 e 2 

294 

490 

15 

P 2 

2 e 1 81e^ 

62e^ 

e 1 37e’^ 

e 1 ISe’^ 

21 245 

525 

35 294 

15 150 

P3 

2 e 229e^ 

58e^ 

37e 227e^ 

e 7e2 

7 735 

525 

105 1470 

15 150 

P4 

8 e lOSe’^ 

0 

lie 229e^ 

e 7e2 

21 735 

35 2450 

15 150 



Thin prolate particles (A —>■ 00 ) 




Total 

Unsteady 

Convective 

Particle 

Pi 

7 

1 

13 

0 

301og2A-45 

81og2A-12 

1201og2A-180 

P 2 

1 

1 

1 

0 

101og2A-15 

81og2A-12 

201og2A-60 

P3 

0 

0 

0 

0 

P4 

0 

0 

0 

0 


(A —?► 0), thin prolate particles (A —>■ oo), and nearly spherical particles. For nearly spherical particles we 
dehne a small parameter e as follows 

A = —^— for prolate spheroids (e > 0), 

1 — e 

A = 1 + e for oblate spheroids (e < 0). 

The asymptotic results for A —>■ 0, A —>■ oo, and for |e| —>• 0 are shown as red dashed lines in Fig. 


B. Linear stability analysis at infinitesimal Re^ 


The effective equations of motion have two special polar angles 9 across which no orbit may pass, 
regardless of the values of /3a. These ang es are 0 = 0 (the vorticity direction) and at 0 = 7 r /2 (the flow-shear 
plane). In the Jeffery dynamics (Re^ = St = 0) the two orbits are called ‘log-rolling’ and ‘tumbling’, and 
they are both marginally stable, just like all other Jeffery orbits. When the /3a are non-zero but infinitesimal, 
the log-rolling and tumbling Jeffery orbits still exist for any finite aspect ratio, but their stabilities change. 

We quantify how particle and fluid inertia lift the degeneracy of the Jeffery orbits by computing the 
stability exponents 7 for the log-rolling (tlr) and tumbling (tt) orbits. The stability exponent is the 
exponential growth rate over one period of the orbit: 


7 = T-i lim log |J 0 (rp)/J 0 o| = T-i 

OC70— 


Alp do 

^89 ’ 


0 


(45) 
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10“^ 10° 10^ 
A 


10 “^ 10 ° 10 ^ 
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FIG. 2. Color online. Coefficients jSa in Eq. ( |42[ ), a = 1..4 as a fnnction of particle aspect ratio A for Re^ = St. 
Solid line shows the snm of all contributions. The other curves show the partial contributions from particle inertia 
(dotted), unsteady fluid inertia (dashed) and convective inertia (dash-dotted). The red dashed lines show the 
asymptotic solutions in Table [II| (first column). Circular markers show result of numerical integration of Eq. (33 1 for 
certain values of A. 


where Tp = 47r/\/l —is the Jeffery period. As Res—i-O we find 

Pi , l-\/l-A2 Pi 

7t = -^H--(A/32-/3i), 7LR = ^ • (46) 

For Res = St these two exponents are shown as function of particle aspect ratio in Fig. Also shown are 
their limiting behaviours in the thin oblate limit (A —>■ 0) 


7 t 1 / 7 34 TttA 

(-53248-b 192007r- 17287r2 - 17287r3 -b5677r4) ;s ^2 


Tlr_ 1_ /5 256 2 

Res ~ 12 ^80 457ry ^12 1357r2 320/ ’ 

in the nearly spherical limit (e —>■ 0) 

7 t 2e 59e^ 7 lr 2e 103e^ 

R^ ^~16«)’ Re7 ^ ~ 2940 ’ 

and in the thin prolate limit (A —>■ 00 ) 

7t _1_ 7lr 1 

R^ 45-30log2A’ R^ Iff^' 


(47) 


(48) 


(49) 


Fig.i shows that prolate spheroids of all aspect ratios are unstable at the log-rolling position, and stable 
at the tumbling orbit. For nearly spherical particles there is a bifurcation: log-rolling and tumbling switch 
stabilities. For oblate spheroids the log-rolling position is stable for any aspect ratio. 
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FIG. 3. Color online. Stability exponents of log-rolling (top panel) and tnmbling (bottom panel) as a function of 
particle aspect ratio A for infinitesimal Res = St. Solid line shows the snm of all contributions. The other curves 
show the partial contributions from particle inertia (dotted), unsteady fluid inertia (dashed) and convective inertia 
(dash-dotted). Red dashed lines show asymptotic results Eqns. (47]|49|. 


For oblate particles there is a second bifurcation at Ac « 1/7.3 where the tumbling orbit becomes stable. 
Clearly, this behavior is caused by the convective inertia of the fluid (see the dash-dotted line in Fig. [^. 
For sufficiently oblate particles both log-rolling and tumbling orbits are stable, and the long-time dynamics 
depend on the initial orientation of the particle. Between the two now stable orbits a new unstable limit 
cycle is born, separating the two basins of attraction. 

Fig. E] shows how the shape of this limit cycle depends upon the particle aspect ratio. Close to the 
bifurcation the limit cycle lies in the neighbourhood of the tumbling orbit. But as A —>■ 0 the limit cycle 
approaches the log-rolling orbit. We have computed the stability exponent of the limit cycle at infinitesimal 
Res by numerically integrating Eqs. (42). The result is shown in Fig. We see that 7 lc > 0, and its 
magnitude is of the same order as that of 7 t. 


V. DISCUSSION 


Effective equation of motion. Eq. (42) is an effective equation of motion for the orientational dynamics of a 
neutrally buoyant spheroid in a simple shear flow. How the dynamics depends upon the particle aspect ratio 
is determined by four coefficients /3i ,..., Pa. Fig. shows the four functions /la(A). Limiting behaviours of 
the Pa are tabulated in Table [H] We see that the /3-coefficients tend to zero as A —)■ oo, but they approach 
constants as A —>■ 0. In both limits the contribution from particle inertia must tend to zero because the 
volume of the particle does. The effects of ffuid inertia vanish as A —>■ oo because the particle effectively 
disappears in the slender-body limit, the perturbation caused by the particle decreases as ~ 1/logA as the 
asymptotic form in Table [H] shows. We remark that the leading-order term in this asymptotic form makes 
a substantial correction to the slender-body theory for aspect ratios of order 30. 

An oblate particle, on the other hand, always presents no-slip boundaries to the fluid, with an area of the 
order of ~ as A —)■ 0. Therefore the contribution of fluid inertia approaches a constant. We note that the 
asymptotic forms of the coefficients Pa listed in Table [n| yield accurate values for A < 1/30 and A > 30, as 
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FIG. 4. Color online. The shape of the limit cycle for different aspect ratios A < Ac = 1/7.3. Trajectories are 
projected onto the unit disk by [X, T] = , n 2 \ (equal area projection). The tumbling orbit is the 

unit circle, log-rolling is the center point. The flow-shear directions are indicated in the background. Parameters 
are, starting from the outer mos t (tumbling) orbit: A = 1/7.2, 1/7.4, 1/8, 1/10, 1/15 and 1/25. Data created by 
numerically integrating Eq. (421 with Res = 10“^. Markers are spaced equally in time. 



FIG. 5. Color online. Stability exp onent quo of the unstable limit cycle as a function of aspect ratio. Gomputed 
by numerically integrating Eqs. (42 1 for Res = 0.05 (solid line). The limit cycle bifurcates at Ac, indicated by an 
arrow and the dotted continuation of the numerical result. 


Fig.i shows. 

We see in Fig. |^that the particle-inertia contribution to the coefficients is always much smaller than 
the fluid-inertia contributions. In general both unsteady and convective fluid inertia contribute, and it would 
be qualitatively wrong to neglect one of these terms. This is due to the fact that the timescale of the particle 
motion is the same as the timescale of the flow, and it raises the question under which circumstances both 
effects may matter for the tumbling of small particles in unsteady flows, and in particular in turbulence. 

Linear stability analysis at infinitesimal Re^. The stability exponents of tumbling and log-rolling orbits 
are shown in Fig. We find that the log-rolling orbit is unstable for prolate spheroids of any aspect ratio, 
tumbling is stable for prolate spheroids, and no other orbit exist at infinitesimal Rej. For moderately oblate 
particles with aspect ratios A > Ac « 1/7.3 the stabilities are reversed: log-rolling is stable, tumbling is 
unstable, and no other periodic orbits exist for infinitesimal Re^. At A = Ac there is a bifurcation where an 
unstable periodic orbit is born close to the tumbling orbit, which in turn becomes stable. As A becomes even 
smaller, the unstable orbit moves closer to the log-rolling orbit (Fig. |^. We remark that the asymptotic 
forms (47) and (49) of the stability exponents yield very accurate approximations for the log-rolling exponent, 
save for aspect ratios close to unity. For the tumbling exponent the asymptotes do not work equally well. 

Our results are in agreement with results of recent DNS studie^SMllSl determining the orientational dy¬ 
namics of a neutrally buoyant spheroid in a simple shear flow. These studies are conducted for a number 
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of different aspect ratios with shear Reynolds numbers ranging from moderate to large. At the smallest 
values of Re^ accessible in the DNS no stable log-rolling is found for prolate spheroids of any aspect ratio. 
For oblate particles with aspect ratio A = 1/5 DNS show stable log-rolling and unstable tumbling at the 
smallest Reg that were simulated^, also in agreement with our results. There are no simulations for particles 
for A < Ac at small Re^. 

SaffmarP predicted that log-rolling is stable for nearly spherical prolate particles, at variance with the 
behaviour described above. We do not know why the original calculation fails to give the correct stability of 
log-rolling. Since no details of the calculation are given it is difficult to figure out the precise origin of this 
discrepancy. Subramanian & KodP^ also computed the stability of the log-rolling orbit for nearly spherical 
particles and came to the same conclusion as Saffman, different from ours. We have compared the small-e 
limit of our calculation to the results of Ref. nzi and find that the particle-inertia correction to the equation 
of motion agrees, Eqs. (3.15) and (3.16) in Ref. [T71 But the fluid-inertia correction does not satisfy the 
symmetries of the problem. We believe that this explains the discrepancy. 

We have independently calculated the stability of log-rolling for nearly spherical particles by expanding 
the particle-angular velocity jointly in e and Re^, using spherical harmonics as a basis selP^. The results 
of this calculation agree to order e with the results presented above. Further we have checked that the 
particle-inertia correction in Eq. (42) is consistent with the results obtained in Ref. |9l We also compared 
the slender-body limit of our resu ts to the prediction of Subramanian & Koch for the dynamics of slender 
fibred and found that the fluid-inertia corrections agree (up to a factor of Stt). 

These observations indicate that the results presented in this paper are correct, explain the results of DNS 
and resolve the puzzle concerning the stability of log-rolling of spheroids in a simple shear at small Re^. 

A new benchmark for DNS at small RCg- Recently a number of groups have developed DNS codes based 
on the lattice Boltzmann method to simulate the dynamics of particles in flow^^lf^. Much effort is spent on 
validating the model, studying for instance the effects changing grid size, time step, size of the simulation 
box, and so forth. The benchmark adopted is often the question whether Jeffery orbits are seen for a 
neutrally buoyant spheroid in a simple shear at small Reynolds numbers. But the limit Re^ = 0 can never 
be strictly reached in the simulations. DNS at small values of Re^ (specifically: in the linear regime), by 
contrast, allow precise comparisons with the results obtained in this paper. One could for instance compare 
trajectories, stability exponents, and period times. We thus expect that our results can serve as benchmarks 
for present and future DNS codes. 


VI. CONCLUSIONS 

In this paper we have derived an effective equation of motion for the orientational dynamics of a neutrally 
buoyant spheroid suspended in a simple shear flow. The equation is valid for arbitrary aspect ratios and to 
linear order in Re^, at small but finite shear Reynolds numbers. The effective equation of motion allows us 
to determine how the degeneracy of the Jeffery orbits is lifted by weak inertial effects. We have determined 
the bifurcations that occur at infinitesimal Re^ as the particle aspect ratio changes. For prolate spheroids 
log-rolling is unstable, for oblate spheroids it is stable. Tumbling in the shear plane is stable for prolate 
particles and unstable for nearly spherical oblate particles. For thin disks with aspect ratios A < 1/7.3, 
both log-rolling and tumbling are stable. An unstable limit cycle separates the basins of attraction of the 
periodic orbits. 

Our results imply that tumbling and log-rolling orbits survive a finite perturbation whose magnitude 
depends on the aspect ratio A. It would be of interest to derive a bifurcation diagram in the A-Reg-plane 
for small Reg. We plan to determine how the small-RCg region of this diagram connects to the intricate 
bifurcation patterns that were found by Rosen, Lundell & AiduiP at larger shear Reynolds numbers. We 
expect that the results summarised here can guide numerical computations with the lattice Boltzmann 
method that become difficult at small Rcg and large aspect ratios. 
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Appendix A: Solutions to Stokes’ equation 


In this Appendix we solve the steady Stokes’ equation for an arbitrarily aligned spheroid in a general 
linear flow u°° = A°^r. The calculation is a special case of the calculation by Jeffery®. However, instead of 
the ellipsoidal harmonics that Jeffery used, we employ a finite multipole expansion, following Chwang and 
Wu^. The purpose of this Appendix is to derive an explicit closed form expression for the Stokes flow field, 
suitable for evaluation in the reciprocal theorem. For a more general description of the method we refer to 
the book by Kim and Karrila^i^. 

Formulation of the problem. Stokes’ equation reads: 

djdjUi = dip , diUi = 0, (Al) 

with no-slip boundary conditions on the surface S of the particle 

Ui = EijkOJjrk for r G S. (A2) 

Here ujj is the angular velocity of the particle. Furthermore it is assumed that the flow remains unperturbed 
at infinitely far away from the particle 


Ui = Ui as r —>■ oo . 


(A3) 


We solve for the disturbance flow m' = Ui — uf° that satisfies Stokes’ equation (Al I with boundary conditions 


4 = Sijki^jrk — u°° for r £ S, u' = 0 as |r| —>■ oo . 


(A4) 


We decompose the linear background flow uf° into its symmetric and antisymmetric parts, defining the 
vector 0“ and strain S°° by 


uT = = Eijk^frk + S°°rp 


(A5) 


Finally, in terms of the ‘slip angular velocity’ Vli — — Wi, the problem to be solved reads 


djdju'i = dip ', diu'i = 0 , 

= -£ijk^jrk - Sf°rj for r £ S, 

u'i = Q as |r|—>- 00 . (A6) 


Multipoles. We solve Eq. (A 61 by a hnite multipole expansiorfiSIlI]. The multipoles are the Green’s function 


for the Stokes’ equation, and its derivatives. In this Appendix we use the shorthand notation Gij^k = dkGij- 
The multipoles needed to solve for the fluid velocity field around particles in a linear flow are 


Sif^Xj Sjj^Xi 3XiXjXj^ 

yij,k — I I 

n QXiXj 

GinM = V -^ , 


G ^(^X'X X 

Gij^llk — ^ Gij,k — c {,^ij^k 


(A7) 
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The following two higher-order multipoles are required in the reciprocal theorem. We include them for 
reference: 


Qij,km — 




^ik^i 




3XjXf-Sim 

3XiX]^Sjm 


SXjXjjiSik , ^Xf-Xm^i 


r 

^X'iX'iYlS ■ 


5 

m^jk 


IbXiXjXkXj] 


Qijjllkm — 


rpO rpO 

‘^OXiXm^jk ^^XiXjSfi^rn 


30XiXf^S jf] 


r' 

30xjXf^Sif] 


30XjXiji^ik 


+ 


r‘ 

SOxj^Xm^i 


'P' p‘ p‘ 

^Sim^jk ^^ik^jm ^^ij^km 


210XiXjXkXm 

J.9 


(A8) 


Note that we use the “Oseen tensor” notation. The Green’s function for the Stokes’ equation is in fact 
Gij = Qij /Stt. It is convenient to split the dipole contribution Qij^k into its antisymmetric (‘rotlet’) and 
symmetric (‘stresslet’) parts. They are 


^ij,k ~ 2 Gikj) 

^ij,k ~ 2 A Gik,j) 


2 i^ik^j j 

Sf^jX'i ^X'iXjXf^ 


(A9) 


Spheroidal multipoles. Whereas the flow around a spherical particle may be represented by multipoles 
anchored at a single po int, r epresenting the flow around a spheroidal particle requires a weighted line 
distribution of multipoleJi^llII. yve therefore dehne the ‘spheroidal multipoles’ as the following distributions, 
note especially the different weights for higher-order multipoles: 

Qij,k{r.ri) = f dC(c2 - ^‘^)GPj,kir “ 

J —C 

Qij,k(r,n) = j dC(c^ - “ ?«), 

= y dC(c2-^n). (AlO) 


The constant c is related to the spheroidal geometry. Prolate and oblate coordinates are obtained by rotating 
an ellipse around its major or minor axis. We call the distance between the foci of the underlying ellipse d, 
and then c = djl for prolate coordinates, and c = id 12 for oblate coordinates (see definition of coordinate 
systems in Appendix [C|) 

In order to write down explicit tensor expressions for the spheroidal multipoles we introduce the integrals 
and AT” by 

2 jn _ Tn+2 

J" - - 2c2C+2 + /"+4. (All) 


r = 


/" = 
= 


The spatial variation of the functions depends upon |rp and r-n only. Further properties and evaluation 
of the integrals are discussed in Appendix]^ With and we express the spheroidal multipoles explicitly, 
for example the spheroidal rotlet: 

— i^ikXj {^ij'f^k ^ik'^j^J^ ■ 




























19 


The integrals play the same part in spheroidal geometry as does l/r"* in spherical geometry. The 
spheroidal stresslet and quadrupole are given by 

^ ~ jTkJ^ ^jk'kliJ^ 

+ ‘i{nirjrk + njTiTk + nkUrj)Jl 
- S{unjnk + TjUiUk + rkriiUj)!^ 

+ 3ninjnkJ^, (A- 12 ) 


= -^{SjkTi + dikVj + 5ijrk)Kl + SOnrjrkK^ 

+ 6{Sjkn^ + S^krij + 5ijnk)Kl 
- 30(rir^nfc + Virkrij + rjrkni)K] 

+ 30{ninjrk + runkr^ + njnkri)K^ 

— SOnirijUkK^ ■ (A-13) 


Solution by a finite multipole expansion. The spheroidal multipoles are functions that satisfy Stokes’ 
equation, and a suitable linear combination of them also satisfies the no-slip boundary condition on the 
surface of a spheroid with symmetry axis n. The remaining problem is to determine the coefficients for this 
linear combination. 

Following Kim and Karrila^ we use the following ansatz for the disturbance flow field; 

'^i ~ Qij^k^jkl [(-^ T ^ i^lm '^1‘^m)) T ^ ^lmn‘^m^no’^o\ 

+ (^Qij,k + oiQ?j,iik) 

X Ujkim B -j- C — C {Sjlmnknm + Sklm'njnm) fli] , 

where 

'^jklm ~ {kljUk — —Sjk){ninm — —Sim) ) 

^jklm ~ kljSkinm “t“ nkSjlHm T njSkmk^l “t“ n.kSjm’kll ^njHkninm : 

^jklm ~ SjkSlm “t” SjlSkm “t” SklSjm 

T SjkniTlm “t“ Slmkljnk njSkinm nkSjinm 

^jSkm'k^l nkSjmkll njnkTlinm ■ 


Given the ambient strain S°°, and angular slip velocity fli = ftj — wf we must determine seven unknown 
scalars, which may depend upon the particle shape: , B^, , and a. When the coefficients 

are known, Eq. (A14) is the sought Stokes solution. _ 

In order to match the linear boundary condition Eq. (A4) we need the combinations of and in the 

ansatz to be constant on the particle surface, much like the scalar function 1 /r™ is in spherical geometry. 
Upon examination, the functions J 3 and are constant on the spheroidal surface. Further, the functions 


J 3 and Kl can be written as J 3 = n,r 


jJ'^i and K\ = 




where and are constant on the 


spheroidal surface. The remaining spheroidal functions and AT" which appear in the ansatz (A14| are 
more complicated. However, it turns out that they appear only in the combinations J" — lOaAT". We 
therefore choose 


lOK^ 


surface 


1 

8 (A 2 - 1 ) ■ 


With this choice of a it holds that, on the surface of the spheroid. 


Jl - lOaAT^ = 0 , 

Jl - lOaAT^ = 0 , 

Jl - A)aKl = Ij 3 - 2a^5^, 

Jl - lOaKl = pi- §aKl , 


(A16) 


(A17) 
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both for prolate and oblate spheroids. 

In order to extract the six independent equations for the six remaining coefficients we exploit that the 
boundary condition must be satisfied for any choice of rij, and S^. First, with = 0 we contract 


Eq. (A4) with n. 


S^rij and finally 


^ 7 ,jj H k and s - 


ipqiiq^pjknj^k- Secondly, with Qj = 0, we contract Eq. (A4) with rii, 
These six equations together have only one solution. We tabulate the 


ij CXJ-J-t-i IIIJ-CIIJ.^ 

resulting expressions for both oblate and prolate spheroids in Table El _ 

Computing the torque on a body due to this flow is straightforward, because by constructiorP^l the torque 
on a body due to the rotlet flow Ui = Q^i.ejki-Ai is = —IfiTrA;, where Ai is the rotlet strength. The 
minus sign is due to the fact that the torque is exerted on body by the flow. To compute the torque from the 
spheroidal rotlet ( A10[ ) we linearly superpose the contributions from all the contained rotlets. The torque 
from the flow Ui = Qijk^jkiBi is therefore 


Ti = -16^ r de(c2 - e)Bi - 

J —C 


647rc^ 


-Bi = C(^Bi. 


(A18) 


The factor cj depends only on the aspect ratio of the particle (see Table III). 


Appendix B: Spheroidal integrals 


In order to solve Stokes’ equation and evaluating the volume integrals in the reciprocal theorem we need 
to solve integrals on the form 


|r|^,r ' 


n)= d^ 


r — ^n\' 


(Bl) 


First, when matching boundary conditions we must evaluate the integrals with r on the surface of the 
spheroidal particle. Second, when evaluating the reciprocal theorem we need to integrate products of two 
or three multiplied with the components of the spatial coordinate r over the entire fluid volume outside 
the particle. Therefore we express the functions in a spheroidal coordinate system with symmetry axis 
x' along n. This is accomplished by a rotational change of variables r' = ilr, x' = i2n, where the latter 
equality defines a rotation R. The absolute value (distance) between r and is preserved by a rotation, 
and the integral is transformed into 


/" = 


dC 


[(a;' - ^2 + (y ^)2 (^7)2]- 


(B2) 


This form is equivalent to the integrals in Chwang and Wu^^. Geometrically, Eq. (Bl| represents a 


fine source along the direction n. The rotation R places the line source along the x'-axis in an auxiliary 
coordinate system. The result is a function of jrp and x' = x' • r' = n • r. 

Explicit expressions for may be found by direct integration, or by a recursion formulsP^. Since we 
require only a finite number of integrals, we simply perform the direct integration once and for all and save 
the result in a table. 

Finally, when evaluating the term corresponding to unsteady fluid inertia in the volume integral of the 
reciprocal theorem, we need to compute the derivatives of with respect to the moving vector n. By 
differentiating Eq. (Bl| we derive the following formula: 


_d_ 

dn. 


— 77 ?V — 77777 ■ 


n+2 


(B3) 


Appendix C: Spheroidal coordinates 

Both oblate and prolate spheroidal coordinates are extensions of a two-dimensional elliptic coordinate 
system (Ci)'? 2 )- The ^i-coordinate represents concentric ellipses, while ^2 represents the corresponding 
hyperbolas. Their intersections give unique coordinates in the x-y-plane. An azimuthal angle of revolution 
(f> denotes the extension into three dimensions. 

Oblate spheroidal eoordinates. Start with the cc-y-plane, and place an ellipse of focal distance d with its 
minor axis along the cc-axis. Now revolve the ellipse by 27r around the a;-axis to produce an oblate spheroid. 
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Then represents concentric oblate spheroidal surfaces, ^2 represents the corresponding hyperbolic surfaces, 
and we call </> the angle of revolution. The coordinate equations are 

X= -66 , 

y = ^\/6 + i\/i-6cos<)i, 

2= + (Cl) 

The coordinate ranges are 0 < 6 < OO) ~1 < 6 ^ 1 0 < ^ < 27r, and the volume element dV = 

(6 + 6) d6d6d<)i. 

In this paper we treat oblate spheroids with dimensionless major axis length unity, and minor axis length 
A. These lengths determine the focal distance d as 


d = 2Vl-A2, 


and the particle surface is parameterised by 


Ap) _ ^ 

61312 - 


(C2) 


(C3) 


Prolate spheroidal coordinates. Start with the x-y-plane, and place an ellipse of focal distance d with its 
major axis along the x-axis. Now revolve the ellipse by 27r around the a;-axis to produce a prolate spheroid. 
Then 6 represents concentric prolate spheroidal surfaces, 6 represents the corresponding hyperbolic sur¬ 
faces, and we call (j) the angle of revolution. The coordinate equations are 

X= -66 , 

y = ^\/^l 

^\/ 6 -l\/^“^ 2 sin(/>, (C4) 

The coordinate ranges are 1 < 6 < OO) < 6 ^ 1 ^nd 0 < ^ < 27r, and the volume element <W = 

(6 - 6) d6d6dti. 

In this paper we treat prolate spheroids with dimensionless major axis length unity, and minor axis length 
1/A. These lengths determine the focal distance d as 


d = 




and the particle surface is parameterised by 




(C5) 


6 a 2 - r 


(C6) 
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TABLE III. Coefficients for Stokes-flow solutions and moments of inertia for prolate and oblate spheroids. These 
coefficients are collected in the book by Kim & KarrilJ^. We tabulate them here for convenience, and because our 
conventions differ slightly from those adopted in Ref. 1101 We remark that some of the coefficients tabulated here 
assume imaginary values. All physical quantities come out to be real-valued. 


Expressions common to both prolate and oblate spheroids 


8(a2-i) 


aR ^ 

4(C-A3-|-A) 




^ Va2-i(a=^-H) 

4(-2CA2-|-C-I-A3-A) 


'. 3/2 


4(-2CA2-|-C-I-A3-a) 




4(2CA2-|-C-3A3-|-3A) 


'. 3/2 


(jS ^ 

2(3C+2A5-7A3-|-5A) 


S _ (a2-i)'’^ [C\+\*-3\^+2) 

~ 8(-2CA2+C+A3-a)(-3CA+A4+a2-2) 


Expressions particular to prolate and oblate spheroids 


C 


Oblate (A < 1) 


Prolate (A > 1) 


-yr^vcot- 






2 \/l - A2 


id 


64 , 


2\3/2 


(1 - A^) 

SttA 


Sa(a^ti) 


2 VA2 - 1 

A 

d 

2 

64^ (A"-l )^^^ 
3A3 


Stt 


15A4 

47r (A^ -h 1) 
l5A4 



























